/*
 * GeoSystem.h
 *
 *  Created on: 22 Aug 2011
 *      Author: Allan
 */

#ifndef GEOSYSTEM_H_
#define GEOSYSTEM_H_

// C++ includes
#include <string>
#include <vector>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include <Core/Multiphase.h>
#include <Core/ReactionSystem.h>

class GeoSystem
{
public:
	GeoSystem();
	
	GeoSystem(
		const Multiphase& multiphase, 
		const ReactionSystem& equilibrium_reactions);
	
	GeoSystem(
		const Multiphase& multiphase,
		const ReactionSystem& equilibrium_reactions,
		const ReactionSystem& kinetic_reactions,
		const vector<string>& primary_species, 
		const vector<string>& equilibrium_species,
		const vector<string>& kinetic_species);
	
	void
	SetTemperature(double T);
	
	void
	SetPressure(double P);
	
	void
	SetComposition(double n);
	
	void
	SetComposition(const VectorXd& n);
	
	void
	SetState(double T, double P);
	
	void
	SetState(double T, double P, double n);
	
	void
	SetState(double T, double P, const VectorXd& n);
	
	const double
	GetTemperature() const;
	
	const double
	GetPressure() const;
	
	const VectorXd&
	GetComposition() const;
	
	const VectorXd&
	GetActivities() const;
	
	const VectorXd&
	GetActivityCoefficients() const;
	
	const VectorXd&
	GetKineticRates() const;
	
	const Multiphase&
	GetMultiphase() const;
	
	const ReactionSystem&
	GetEquilibriumReactions() const;
	
	const ReactionSystem&
	GetKineticReactions() const;
	
	const vector<string>&
	GetPrimarySpecies() const;
	
	const vector<string>&
	GetEquilibriumSpecies() const;
	
	const vector<string>&
	GetKineticSpecies() const;
	
	const double
	Moles(const string& species) const;
	
	const double
	Mass(const string& species) const;
	
	const double
	Activity(const string& species) const;
	
	const double
	ActivityCoefficient(const string& species) const;
	
	const double
	KineticRate(const string& species) const;	
	
	const double
	Acidity() const;
	
	const double
	IonicStrengthEffective() const;
	
	const double
	IonicStrengthStoichiometric() const;
	
	const double
	PartialPressureGas(const string& gas) const;
	
private:
	/// The multiphase object of the geochemical system
	Multiphase multiphase;
	
	/// The equilibrium-controlled reactions of the geochemical system
	ReactionSystem equilibrium_reactions;
	
	/// The kinetically-controlled reactions of the geochemical system
	ReactionSystem kinetic_reactions;
	
	/// The primary, equilibrium and kinetic species of the geochemical system
	vector<string> primary_species, equilibrium_species, kinetic_species;
	
	/// The temperature and pressure of the geochemical system
	double T, P;
	
	/// The composition (moles) of the species of the geochemical system
	VectorXd n;
	
	/// The activities of the species of the geochemical system
	VectorXd a;
	
	/// The activity coefficients of the species of the geochemical system
	VectorXd g;
	
	/// The kinetic rates of comsumption/production of the species of the geochemical system
	VectorXd r;
	
	/// The stoichiometric matrix containing the stoichiometry of every species in the kinetic reactions
	MatrixXd uki;
	
	/// The molecular weight of all species in the multiphase instance
	vector<double> molecular_weights;
	
private:
	void
	InitializeStoichiometricMatrixUki();
	
	void
	InitializeMolecularWeights();
};

std::ostream&
operator<<(std::ostream& out, const GeoSystem& geosystem);

#endif /* GEOSYSTEM_H_ */
